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A Virtual Element Method for elastic and inelastic 
problems on polytope meshes 

L. Beirao da Veiga * * C. Lovadina } D. Mora 


Abstract 

We present a Virtual Element Method (VEM) for possibly nonlinear elastic and 
inelastic problems, mainly focusing on a small deformation regime. The numerical 
scheme is based on a low-order approximation of the displacement field, as well as 
a suitable treatment of the displacement gradient. The proposed method allows 
for general polygonal and polyhedral meshes, it is efficient in terms of number of 
applications of the constitutive law, and it can make use of any standard black-box 
constitutive law algorithm. Some theoretical results have been developed for the 
elastic case. Several numerical results within the 2D setting are presented, and a 
brief discussion on the extension to large deformation problems is included. 


1 Introduction 

The Virtual Element Method (VEM), introduced in [2], is a recent generalization of the 
Finite Element Method which is characterized by the capability of dealing with very 
general polygonal/polyhedral meshes and the possibility to easily implement highly 
regular discrete spaces HUE!. Indeed, by avoiding the explicit construction of the 
local basis functions, the VEM can easily handle general polygons/polyhedrons without 
complex integrations on the element (see [3] for details on the coding aspects of the 
method). The interest in numerical methods that can make use of general polytopal 
meshes has recently undergone a significant growth in the mathematical and engineering 
literature. Among the large number of papers, we cite as a minimal sample ma 01231 
mil [T21 1501 EH [T7] [TO]. Indeed, polytopal meshes can be very useful for a wide 

range of reasons, including meshing of the domain (such as cracks) and data (such as 
inclusions) features, automatic use of hanging nodes, use of moving meshes, adaptivity. 

In the framework of Structural Mechanics, recent applications of Polygonal Finite 
Element Methods, which is a different technology employing direct integration of com¬ 
plex non-polynomial functions, have shed light on some very interesting advantages of 
using general polygons to mesh the computational domain. This include, for instance, 
the greater robustness to mesh distortion H3], a reduced mesh sensitivity of solutions 
in topology optimization H3 ED], better handling of contact problems [7] and crack 
propagation [23]. Unfortunately, Polygonal Finite Elements suffer from some serious 
drawbacks, such as the strong difficulties in the three dimensional case (polyhedrons) 
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and in the use of non convex elements. On the contrary, the VEM is free from the 
above-mentioned troubles, and thus it represents a very promising approach for Com¬ 
putational Structural Mechanics problems. 

Aim of the present paper is to initiate the investigation on the VEM when applied 
to non-linear elastic and inelastic problems in small deformations. More precisely, 
we mainly focus on the following cases: 1) non-linear elastic constitutive laws in a 
small deformation regime which, however, pertain to stable materials; 2) inelastic con¬ 
stitutive laws in a small deformation regime as they arise, for instance, in classical 
plasticity problems. We remark that we are not going to consider here situations with 
internal constraints, such as incompressibility, which require additional peculiar nu¬ 
merical treatment. Virtual elements for the linear elasticity problem where introduced 
in 0Ej. The scheme in the present paper is one of the very first developements of 
the VEM technology for nonlinear problems, and it is structured in such a way that 
a general non linear constitutive law can be automatically included. Indeed, on ev¬ 
ery element of the mesh the constitutive law needs only to be applied once (similarly 
to what happens in one-point Gauss quadrature scheme) and the constitutive law al¬ 
gorithm can be independently embedded as a self-standing black-box, as in common 
engineering FEM schemes. Therefore, in addition to the advantage of handling general 
polygons/polyhedra, the present method is computationally efficient, in the sense that 
the constitutive law needs to be applied only once per element at every iteration step. 
The risk of ensuing hourglass modes is avoided by using an evolution of the standard 
VEM stabilization procedure used in linear problems. However, we highlight that the 
proposed method is described for general d-dimensional problems {d = 2,3), but the 
performed numerical experiments are confined to the two dimensional setting. 

A brief outline of the paper is as follows. In Section [2] we describe the continuous 
problems we are interested in. In particular, we distinguish between the elastic, possibly 
non-linear, case (Section 2.1), and the general inelastic case (Section 2.2). Section [3] 


deals with the VEM discretization. After having introduced the approximation spaces 
and the necessary projection operators (Section 3.1), we detail the discrete problems 
for the elastic case in Section [3.2| and for the inelastic case in Section 3.3 In Section |4j 
combining ideas and techniques from [14] and [2j, we provide some theoretical results 
concerning the convergence of the proposed scheme in the elastic situation. We remark 
that our analysis is confined to cases where the non-linear costitutive law fulfills suitable 
continuity and stability properties, as stated at the beginning of the Section. Section [5] 
presents several numerical examples which asses the actual behaviour of the proposed 


scheme. In Sections 5.1 and 5.2 we consider non-linear elastic cases, while in Section 5.3 
a von Mises plasticity problem with hardening is detailed. Furthermore, an initial 
brief discussion about a possible extension to large deformation problem is included 
(Section [5~4| ). Finally, we draw some conclusion in Section |6j 

Throughout the paper, we will make use of standard notations regarding Sobolev 
spaces, norms and seminorms, see [8], for instance. In addition, C will denote a constant 
independent of the meshsize, not necessarily the same at each occurrence. Finally, given 
two real quantities a and b, we will write a < b to mean that there exists C such that 
a < Cb. 


2 The continuous problems 

In the present section we describe the problem considered in this paper. Although the 
elastic case could be considered as a particular instance of the inelastic case, we prefer 
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to keep the presentation of the two problems separate. This will allow us a clearer 
presentation of the ideas of the virtual element scheme in the following section. 

2.1 The elastic case 

We consider an elastic body £l cM. d (d = 2,3) clamped on part T of the boundary and 
subjected to a body load f . We are interested, assuming a regime of small deformations, 
in finding the displacement u : Q —>• M. d of the deformed body. 

We are given a constitutive law for the material at every point x € £1, relating 
strains to stresses cr, through the function 

<T = a(x, Vu(®)) € Rsymm C 1 ) 

where Vu represents the gradient of the displacement u. 

Given the law 0> the deformation problem reads 

{ — div a = f in 12, 

u = 0 on T, (2) 

crn = 0 on <912/T, 

where n denotes the unit outward normal to d£l. 

Let now V denote the space of admissible displacements and W the space of its 
variations; both spaces will, in particular, satisfy the homogeneous Dirichlet boundary 
condition on T. The variational formulation of the elastic deformation problem reads 

{ Find u£ V such that 

f cr(x, Vu(x)) : Vv(x) dx = j f(x)-v(.x)dx Vv E W. 

Jo. -hi 

Remark 2.1. The generalization of the results of the present paper to other type of 
loadings (for instance in the presence of boundary forces) and boundary conditions (for 
instance in the presence of enforced displacements) is trivial. Our choice in allows 
to keep the exposition shorter. 

2.2 The inelastic case 

We assume a small deformation regime and restrict ourselves to rate independent in¬ 
elasticity. We consider a material body 12 C (d = 2,3) clamped on part T of the 
boundary and subjected to a body load f(f, x) depending also on a pseudo-time variable 
t € [0,T]. The interested reader can find more details in [231 [27] , for instance. We are 
interested in finding the displacement u : 12 —>• of the deformed body at a given final 
time T. 

We are given an inelastic constitutive law for the material, relating strains to stresses 
a, through the function 


a = a(x, Vu (x),H x ) <E Mf y x ^ m (4) 

where the vector H x contains all history variables at the point x. 

The above rule is to be coupled with an evolution law £ for the history variables in 
time 

Tix = C(x, Vu(i), Vu(x),H x ), (5) 
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where, as usual, a dot above a function stands for a pseudo-time derivative. Since we 
consider a quasi-static problem, at each time instant the stresses and displacements 
must satisfy the equilibrium and boundary conditions in 

We here avoid to write a rigorous variational formulation for the problem above, 
and limit ourselves to the minimal setting that will be needed to introduce the asso¬ 
ciated discrete problem. As in the elastic case, let V denote the space of admissible 
displacements and W the space of its variations. Then, assuming an initial value for 
the history variables, the quasi-static inelastic deformation problem can be written as 

f For all t G [0, T\ find u (t, •) G V such that 

| [ a(x, Vu(t, x), H x (t)) : Vv(x) dx = [ f(f, x) ■ v(x) dx Vv G W, 

where the displacements and history variables are sufficiently regular in time and must 
satisfy the evolution law © almost everywhere. 


3 The virtual element approximation 

In the present section we introduce the virtual element discretization of problems ^ 
and (|6]). In what follows, given any subset lo of (d = 2,3) and k G N , we de¬ 
note by Vy-{u) (respectively 'Pjfc(w)) the scalar (respectively vector with d components) 
polynomials of degree up to k on oj. 

3.1 The virtual spaces and operators 

We consider a mesh P^ for the domain P, made of general polygonal/polyhedral con¬ 
forming elements. For the time being, we only assume that such mesh is compatible 
with the boundary conditions, i.e. that T is union of faces (edges) of the mesh. We 
denote by E G P^ the generic element of the mesh and by / the generic face (or edge if 
d = 2). The symbols He and \E\ will represent, respectively, diameter and volume (or 
area) of the element E. As usual, h will indicate the maximum element size. 

We start by introducing the discrete virtual space for displacements, that is essen¬ 
tially the same as in [3j. We first consider the two dimensional case. Given any E G P^, 
let the local virtual space 

V hlE ■■= {v G [H\E) n C°(E)] 2 : Av = 0 in E, v|, G Vi (/) V/ G dE}, (7) 

where A denotes the component-wise Laplace operator. The space V^e is a space of 
harmonic functions that on the boundary of the element are piecewise linear (edge by 
edge) and continuous. Such space is virtual in the sense that is well defined but not 
known explicitly inside the element. 

Note that 'Pi(E) C Vh t E] in the case of a triangular element, we recover exactly the 
standard Vi space. It is easy to check [3] that a set of degrees of freedom for the space 
Vh,E is simply given by the collection of the vertex values: 

• Pointwise values {v^)}^^ with v denoting a vertex of E. 

Once the above degrees of freedom values are given, since v G Vh.E is linear on each 
edge, the value of v on the boundary dE is completely determined. Therefore, an 
integration by parts allows to compute the integral average of the gradient 

t4t / Vv dx=^-^ /v®n / ds VvGVft >E , (8) 
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with n/ indicating the outward unit normal at each edge /. 

We now define the virtual local spaces for the three dimensional case. Given a 
polyhedron E £ J4, any face / £ dE is now a polygon. We denote by Vhj the virtual 
bi-dimensional space 0 on the polygon / adjusted with three components: 

Vh,f ■■= {v £ n C°(/)] 3 : Av = 0 in /, v| e £ Vi(e) Ve £ df}, (9) 

where the symbol e represents the generic edge of the polyhedron and A denotes the 
planar laplacian on /. We then define 

V h>E = {v £ [H\E)f : Av = 0 in E, v\ f £ V hJ Vf £ dE}. (10) 

The space Vh,E is a space of harmonic functions that on the boundary of the element 
are continuous and, on each face, functions of 14,/- Note that, as a consequence, the 
functions of Vh.E are linear on each edge of the polyhedron. 

Again we note that 'Pi(E) C V^e] i n the case of a tetrahedral element, we recover 
exactly the standard Vi space. It is easy to check that a set of degrees of freedom for 
the space V^e is again given by 

• Pointwise values {'v(^)}u£dE with v denoting a vertex of E. 

An integration by parts exactly as in ([8]) allows to compute, for all E £ 14 the integral 
average of the gradient, provided one is able to compute the face integrals fj v <g) rif ds 
for all / £ dE and v £ Vh.E- Such face integrals can be easily computed by introducing 
the virtual space modification proposed in [T], that we do not detail here. The result is 

/ v<8> n/ = oj u v(v), 
v£dE 

where the scalars {ui l ,} u &dE are the weights of any integration rule on the face that is 
exact for linear functions. 

Once the local virtual spaces are defined, all that follows holds identically in two 
and three dimensions. We can now present the global virtual space 

14 := {v £ V : v| E £ V h , E VE £ 14}- 

A set of degrees of freedom for Vh is given by all pointwise values of v on all vertices 
of £lh, excluding the vertices on T (where the value vanishes). 

In the following, we will denote by II 0 the tensor valued L 2 projection operator 
on the space of piecewise constant functions and by 11^ its restriction to the generic 
element E £ fMore precisely, for any G £ ( L 2 (£l)) dxd , we have (II 0 G)£ = II^(G|e) 
with the local operators defined as 

U °e G \e = / G VEeQ h . (11) 

l-E'l JE 

We have the following important remark, which is a direct consequence of (|8|. 

Remark 3.1. For all functions v £ Vue and. all elements E £ 14, the operators 
n ° E (Vv ) are explicitly computable. 
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We moreover introduce a second projection operator II V , defined on 14 as follows. 
For any v G 74, we have (II v v)£ = n^vj^) G "Pi (E) with the local operators defined 
as 


v(nI(v| s )) = n°;(Vv| B ), 

Y ( u e v ){^) = Y 

u^dE v^lOE 


for all E in fNote that, by definition, II v v is a (discontinuous) piecewise linear 
function on I7/i. On each element E, 11^(v|^;) is the unique linear function such that: 


1. its (constant) gradient equals the mean value over E of the function Vv; 

2. its vertex value average equals the vertex value average of v. 


We notice that the second condition in (12) is only to fix the constant part of II v v 
on each element. Recalling Remark 


3.1 


it is immediate to check that the operator II V 


is explicitly computable. 


3.2 The elastic case 

The main missing step is to introduce the local forms that will be used in the dis¬ 
crete variational formulation. We assume that the constitutive law ([!]) is piecewise 
constant with respect to the mesh 17^. Therefore, instead of er(x, Vu(x)), we will write 
<x£:(Vu(x)) to represent the constitutive law on E, E G and x G E. In addition, foe 
every pair v G V and w G W, we introduce the bilinear forms a E (y, w) and a(v, w) as: 


Therefore, it holds 


ag(v,w) = / <r(x, Vv(r)) : Vw(i) dx, 
Je 

a(v, w) = / <x^(Vv(x)) : Vw(.x) dx. 

J n 


a(v,w)= Y «e(v,w) 


(13) 


(14) 


and, recalling (§, the elastic problem can be written as 

Find u G V such that 

a(u, v) = / f(x) • v(x) dx Vv G W. 
Jn 


(15) 


We now consider, for all E G 17^ and all v/^w/j G 14,Ej the following preliminary 
form 


ah,E(vh, Wft) = [ a E {Il 0 E {Vv h )(x)) : (U° E (\/w h )(x)) dx 
Je 

= \E\ tr E (n° E (Vv h )) : n»(Vw A ,), 


(16) 


where the identity above follows since all the involved functions are constant on the 
element. The above form is Pi-consistent, in the sense that it recovers exactly the 
original form whenever the first entry is a linear polynomial. Indeed, it follows from ( pTTj ) 
and (12) that 
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ahMV’Vh) = [ cr B (II^(Vq)(a:)) : (U° E (Vw h )(x)) dx 
Je 

= [ o\e(V q(x)) : (n° E (Vv/ h )(x)) dx = [ <x E (Vq(x)) : Vw h (x) dx ( 17 ) 
Je Je 

= a E { q,v ft ) Vq.eVi(E),Vv h eV htE - 

However, unless the elements are triangular/tetrahedral, the form c4,e(-, ■) has a 
non-physical kernel that may lead to spurious modes in the solution. We therefore 
follow the idea proposed initially in [2] and introduce the discrete bilinear form 


Sh,E ■ Vh,E x Vh,E —t R, 

Sh,E(y h ,vfh) = h d E 2 ^2 v h (v)vr h (v) Vv h , w h G V h)E - 

uedE 

As discussed in under suitable mesh regularity assumptions detailed in Section 
[4j there exist positive constants c*, c* independent of the element such that 

c* f ||V sym v/ l || 2 dx < S htE (vh,Vh) <c* f ||V sym v h || 2 d* (19) 

Je Je 

for all v/j G Vh,E with n^v/! = 0. In other words, on the orthogonal complement of 
Vi(E) with respect to V^e, the bilinear form Sh,E (■) behaves as the local energy of 
a linearly elastic body with unitary material constants and is thus suitable to stabilize 
ah, E {-, •) form in such case. In order to take into account different material constants 
and also nonlinear materials, the form S^ E {-,-) needs to be multiplied by a positive 
constant a E that may depend on the discrete solution. 

We therefore introduce the following local virtual form on V^e- For all E G and 
all s h ,v h ,w h e V h)E 

ah,E( s hl v h, w ft ) = a h>E {\hi w h) + aE{s h )S htE (\r h - Uj\r h ,w h - U^w h ), (20) 


where the stabilizing parameter a E > 0 depends on the additional entry s/, . We remark 
that the bilinear form ah, E {-\ •, •) is still "Pi-consistent. This follows from © and the 
observation that q — n^q = 0 for every q G “Pi- The choice that we here propose for 
the parameter a E is 

a E (s h ) = |||^(n^Vs h | B )||| ME G Q h , Vs h G V h , (21) 


representing any norm on the fourth order tensor space, for instance the 


with 11 

maximum of the absolute values of all the entries, see Remark 3.2 
We present also the global form 


ah(s h \ v h , w ft ) = a h, E (s h ;v h ,'w h ) Vs h ,v h ,w h G V h . 

Given s/j G a possible virtual discretization of Problem © is 

{ Find u/j G 14 such that 
ah{s h ] u h ,v h ) =< f, v h > h Vv h G V h . 

Above, the load approximation term 

<f ,v h >h= ^2 u u f(v)v(is) 

v&dE 


( 22 ) 


( 23 ) 
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is a vertex-based quadrature rule with weights chosen to provide the exact integral 
on E when applied to linear functions. Furthermore, a reasonable choice for s/j could 
be s h = u h . 

We instead propose a modification of (23), that is more practical from the im¬ 


plementation viewpoint. We assume the usual incremental loading procedure for the 
solution of the nonlinear discrete problem: given a positive integer N, let the partial 
loadings f” = (n/N )f for all n = 1,2,..., TV. Then, given the initial displacement 
(for instance the zero function), one applies for n = 1,2,...,N the iterative procedure 


Find u£ E V) t such that 
^(u^ju^Vft) =< f n ,v h > h Vv h € V h . 


(24) 


The final solution is then u/, = u^. The nonlinear problems above can be solved with 


the Newton scheme. Note that, since the stability constants ag (see (201) are computed 
by using u)) -1 , the tangent matrix in the Newton iterations turns out to be simpler. 
Since N is typically taken large (at least 10, but often much more) the effect of such 
modification is not detrimental for the discrete approximation; the constants cxe are 
only used as scaling parameters and do not enter the accuracy of the algorithm. 

We close the section with some observations regarding the local forms used in 
the scheme. First, we recall that the proposed forms are 7 7 i(i?)-consistent, in the sense 
that for all E £ klh, we have: 

aVE^q, v h ) = [ (Tg(Vq) : Vv^ dx Vs^, v h £ V h E , Vq E V\(E). (25) 

Je 


Identity (25) is a fundamental condition for approximation and, in particular, guaran¬ 


tees the satisfaction of the patch test. Moreover, such forms are explicitly computable 
for any polygonal/polyhedral element (even non-convex). Finally, the constitutive law 
needs to be computed only once per element and thus the method, from this point of 
view, is as cheap as finite elements with one point gauss integration rule. This obser¬ 
vation has an even bigger impact in the inelastic case, where the constitutive laws are 
typically more expensive to compute. 


Remark 3.2. The motivation for choice (21) and (24) is to better mimic the sta¬ 


bility properties of the material for the current displacement. For materials in which 
the stress-strain incremental relation does not depend too strongly on the value of the 


current displacement, the constants cue can be taken as independent of u/ . For in¬ 
stance, a scaling directly proportional to the local material constants could be used. On 


the other hand, the choice proposed in (21) and (24) give good results for a wider range 


of materials. Examples and investigations in this direction can be found in Section [5| 


3.3 The inelastic case 

We start by introducing a sub-division of the “time” interval [0, T] into smaller intervals 
[t n -i,t n ] for n = 1,2, ...,N, where for simplicity we assume that t n = nT/N. We will 
denote the partial loadings by f n = (n/N )f for all n = 1, 2,..., N. 

We assume, as in standard engineering procedures, a constitutive algorithm that 
is an approximation of the constitutive and evolution laws Q, (|5]). In Finite Element 
analysis, this pointwise algorithm can be coded independently from the global FE 
construction and can be regarded as a “black-box” procedure that is applied at every 
Gauss point and at every iteration step. In the present Virtual Element method, we 









want to keep the same approach; in other words, our scheme will be compatible with 
any black-box constitutive algorithm that follows in the general setting below and that 
can be imported from other independent sources. 

We assume that the constitutive law is piecewise constant with respect to the mesh 
f2/!. Let <r e represent the constitutive algorithm for the element E £ 12/j. For any 
x £ E, given a value for the displacement gradient Vu^ -1 (x) at time f n _i, a value 
H x ~ l for the history variables at time t n - \ and a tentative value for the displacement 
gradient Vu^(x) at time t n , the algorithm computes the stresses (and updates the 
history variables) at time t n . We thus write the computed stress as 

<r E (Vu 

As part of the approximation procedure of our method, we assume that the history 
variables H x are piecewise constant with respect to the mesh, and therefore write H E 
to represent the value assumed on the element E £ H/i at time t n . Consistently, H n 
will represent the collection of all {H E }EeCi h - 

In our scheme, instead of applying the constitutive algorithm at Gauss points, we 
make use of the projections introduced in the previous sections and of the same stabi¬ 
lization as in the elastic case. The Virtual Element scheme reads, for n = 1,2,..., N: 

( Find U/J £ Vh (and the updated H n ) such that 

1 a h (u%- 1 ,v%,H n ~ 1 ,v h ) =< f n , v ^ > h Vv ft € V h , 

where the form 


«/»« \v h )= J 2 a hM u h 1 ^ u h^ H E 1 ^h) 

with, for all E E 

a h , E {v n h -\<,H n E -\v h ) =\E\ a E (U° E Vu n h 

+ oi E {nl~ x )s hjE {K -n^ug, Vft -n^v ft ). 

Here above, the bilinear form Sh, E and the scalar a E are calculated as already shown 
in (18) and (21), respectively. Note that, as already mentioned in Section 3.2 the 
constitutive algorithm needs to be applied only once per element. 


h u e ~\ n u B V<) : n u E (Vv,) 


t° 

l E 

T V, 


4 Theoretical results 

We here develop an error analysis for the method described in Section [3.2[ under some 
additional hypotheses on the function cr(x, Vu(x)) = <te(Vu(x)). More precisely, we 
assume that the following properties are satisfied. 

Hypotheses (RPC) 

• The function r >-)■ cr E (r) belongs to C 1 (M dxci ) for every E £ 

• for every E £ the differential (r) satisfies 

1. there exists C a > 0 such that 

^(t) s : s > Callsl) 2 VsGM dxd , (27) 

OT 
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2. there exists Cm > 0 such that 


da E 

8t 


(t) s : t < Cm \ |s| 


Vs,t e 


bdxd 


(28) 


We moreover explicit here the shape regularity conditions that are needed for the 
theoretical results of the present paper. We assume that there exists a positive constant 
C s such that all the elements E of the mesh sequence are star shaped with respect to 
a ball of radius p > C s h E and that all the edges e of E have length h e > C s h E . 

Lemma 4.1. Let the bilinear forms a E (-,-), ah, E (-',-,-) and, •, -) be defined 

by (13), (20) and (22). Suppose that the Hypotheses (RPC) introduced above are satisi- 
fied. Then, it holds 


|v/j - n < a h (s h -, v h , v h - w h ) - a h {s h - w h , v h - w ft ) Vv/,, w ft , s h € V h . (29) 

a E (v,r) - oe(w, r) < |v - w| lijE |r|i ijE Vv,w,reV. (30) 

a,h, E (sh\ v h ,r h ) - a h , E (s h \ w h ,r h ) < \v h - Wh\i, E \rh\i, E V v h , w h , s h , r h E V h . (31) 


Proof. We first note that (|27|) and (|28|) , together with (|21|), imply the existence of 

(32) 


positive constants ci and C 2 such that 
ci < a E (s h ) < c 2 


E n h , V S/) . E v h . 


Step (i): proof of (29). From (27), we deduce that 

(cr E ( s) - cr E (t)) : (s - t) > C||s - t || 2 
Therefore, for every v, w E V we have 


Vs,t E 


bdxd 


(33) 


a E (v, v — w) — a E (w, v — w) = / (cr E (Vv) — <te(Vw)) : (Vv — Vw) > C|v — w| 2 E , 

Je 

(34) 

by which 


v — w 


l,f2 


< 


a(v, v — w) — a(w, v — w) Vv,w£V. 


(35) 


For every v h ,w h ,s h E V h , we have (see fl20|)) 

ah,E(sh 5 Vh, Vh - Wh) - a h)E (s h ; w h , v h - w ft ) 

= a h , E (v h , Vh - Wh) - a h , E { Wh, v h - w ft ) 

+ a E (s h )S hjE ((v h - w/ t ) - n|(vh - w/ t ), (vh - Wh) - n^(vh - w,,)). 


(36) 


We now notice that (see (16)) 

5 h,E(vh, Vh - wh) - a htE ( w /M Vh - Wh) 

= /' <T E (n° B (Vvh)): (n^(Vvh) - n° E (Vwh)) 

Je 

- f a E (U 0 E (V Wh)) : (n^(Vvh) - n^(Vwh)) 

J E 

= I [o- E (n^(Vvh)) - <r E (4(Vwh))]: (n^(Vvh) -n^(Vwh)). 

J E 

First using (33) with s = Ilg(Vvh) and t = n^(Vwh), then recalling we get 
ah,e(vh,v/i - Wh) ~a hiE (w h ,v h - w h ) > C'||n|(Vvh) - n^(Vwh)||o ;E 

= c || v(nl(vh - W h))lliU = c 1 nlK - Wh)|? iE . 


(37) 


(38) 
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( 39 ) 


In addition, we have, using fl32ft and (19): 


aE(Sh)Sh,E(( V h ~ W ft) ~ n I( V ^ “ W h), (Vh. ~ W h ) - II %(v h - W/J) 
> C\(y h - w h ) - Il|(v ft - w h )\j E 


Combining (36) with (38) and (|39j) , we infer 

\vh-Wh\l,E ^ a ht E(s h ;v h ,v h -w h )-a ht E(sh;w h ,Vh-w h ) 


Summing up over all the elements, we get 
I Vh - Wh\i,n < a h (s h ;v h ,v h - w h ) - a h (s h : w h , w h - w h ) 


Vv h ,w h ,s h G V h . (40) 


Vv/j, w/j, Sft G V/j. (41) 


(**); proof of (30) and (31). From ( |28| ), we deduce that 

(cte(s) — CT£(t)) : t < C||s — 111 ||r|| Vs,t,r G R dxd , 


from which we easily get (30): 


ob(v, r) - oe(w, r) < |v - w|i )E |r| 1)E 


V v, w, r G V. 


(42) 


(43) 


We now notice that (see ©) 

dh,E(v h ,r h ) - d h)E (y*h,*h) = [ [<r E (II^(Vv^)) - o- E (II^(Vw^))] : U° E (Vr h ). (44) 

J E 


Using (|42j), identity (|44j) yields 

dh,E(,v h , r h ) -ah, E (wh,r h ) < |v fe - w ft |i )E |r fe |i )B \/v h ,w h ,r h G V h . 


(45) 


To continue, since Sh,E(', ■) is a bilinear form and using continuity arguments, we have 
for every s/,, G 14 (see 

atE(sh)S htE (vh ~ n|(v fe ), r h - n|(r ft )) - a E (s h ) S hiE (vr h ~ V^(w h ),r h - Il|(r h )) 

= OLE(s h )S h , E ((v h - W h) - U^(\r h - w h ),r h - Il|(r fe )) 


< 


V/, - Wh\l,E\ r h\l,E- 


From (20), using (45) and (46), we deduce (31). 


(46) 

□ 


Theorem 4.1. Let u G V be the solution of Problem (§• Given any s h G 14, let 
u h G Vh be the solution of Problem (23): 


{ Find u h G 14 such that 
a>h{s h ',u h ,v h ) =<f,v h > h Vv/j G 14- 

For any u/ G 14 and u T G L 2 (U) such that u^\ E G Vi(E), it holds: 


(47) 


< f,V h > h -(f,Vfc) |,| | /.ON 

lu-Ufcl^Q^ sup -j—j- b |u - u/^n + |u - u^i^, (48) 

Vhe\4 l v h|i,n 


where (•,•) denotes the [L(Ll)] d -scalar product. 
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Proof. Given u/ E 14, we set 5 / l = — u j. For every E L 2 (S7) such that u^\ E E 

Vi (E), using (29) we have 


|u h ~ u/l^n < ah(s h ; u h , S h ) - a h (s h ; u/, d h ) 

=<f,5h>h~ ^2 a h,E(Sh]U-I,8h) 

E&Q h (49) 

^ f j 5ft P h ^ [P'hjE^h^ U/> 5ft) ^h,Eiphi j 5ft)] 4“ ^h^iph 1 ^7r j 5ft) ^ • 

E(zQh 

Since (|25]) implies a h ,E(Ph\ u f , S h ) = a E (u n ,S h ), from (|49]> we get 

|Uft - u /|l,a ~ <f ,$h>h- J2 { [ a h,E{ s h\ u /j 5ft) a hjE (s h \ IIttj 5ft)] + a E (u n , Sh) ^ 

^ f 5 >h ^ ^ {p'hjE^hi 5 &h) ^h,E{^hi ^7r? ^/i)] 

EtzClh 

- J2 [ a E(n n ,S h ) - a E {u,S h )] - a(u,S h ) 

E&n h 

= [ < f, 5ft >h —(f, Sh)] — ^2 [ a h,E^>h\ U/, Sh) — a h ,E( Sh', Utt, 5ft,)] 

Etzflh 

- 5Z [aE(vL n ,S h ) - ci E (u,S h )}. 

E€.£lh 

(50) 


We then obtain, using (30) and (31) 


I ,2 / ( < f,Vft >ft -(f,Vft) | | it / K1 v 

I u ft-u/| in < sup -;—;-(- |u/ - u^-lm + |u*. - u|i : n I |dft|i,n, (51) 

\ v ftGVft |Vft|l,n J 

by which, recalling that 5ft = Uft — u/, we infer 


Uft-u/|i,n< sup 

VftGVft 


< f, Vfe >h ~(f,V/t) 

|v/i|i,n 


+ |u/ - UttIpo + Iutt - ulpn. 


(52) 


The triangle inequality thus gives 


< f,Vft >ft (f, Vft) III I /roN 

|u-Uft|l,n< sup -;—j-(- |u - U/^fi + |u - U^i^. (53) 

VftGVft | v h|l,f2 

□ 


Remark 4.1. Theorem f.l applies also to Problem (24) at the final step N. Indeed, 
it is sufficient to make the choices f = f iV , s at = u /( _ in (47) , and to identify lift 


in (|47j) with in ([24]). 


Corollary 4.1. Following the same notation of Theorem f.l. let moreover u E [H‘ 2 (H)] d . 
Then the linear convergence bound holds 


u - Uftlpn < h |u| 2 ,q. 


Proof. The results follows immediately combining Theorem |4.1| with standard polygonal 
approximation estimates for the spaces V^e an( l 'Pi(-E), see [2, 125] , □ 
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Figure 1: Sample meshes: £1^ (left) and (right). 



Figure 2: Sample meshes: (left) and (right). 


5 Numerical tests 


In the present section we test our virtual method. In the first two examples (see 
Sections 5.1 and 5.2), the body occupies the region fl := (0, l) 2 , where lengths are 
expressed in meters. We employ the following types of mesh (see also Figures [l]j2]): 


• Structured hexagonal meshes. 

• fl 2 : Non-structured hexagonal meshes made of convex hexagons. 

• : Regular subdivisions of the domain in IV x IV subsquares. 

• Trapezoidal meshes which consist of partitions of the domain into congruent 
trapezoids, all similar to the trapezoid with vertexes (0,0), (^,0), (^,|), and 

(o. s)- 


In what follows, Nh denotes the number of vertices in the mesh under consideration. 
To test the convergence properties of the methods, we introduce the following dis¬ 
crete maximum norm: for any sufficiently regular function v, 


|||v|||o,oo := max |v(v)|oo 

vEVh 

where Vh represents the set of vertexes of Qh and | • |oo denotes the l°° 
We also introduce the following discrete H 1 like norm: 


v 1,2 


■~ \ ^2 h e 




dv 

die 


1/2 


0,e , 


(54) 

vector norm. 

(55) 
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where £h and h e denote the set of edges in the mesh and the length of the edge e, 
respectively. Moreover, t e denotes one of the two tangent vectors to the edge e, chosen 
once and for all. Accordingly, we denote by 

E 0 ,oo '■= lllu-uhllloo £i ) 2 := |||u-u fc |||i i2 

the corresponding errors and we measure the experimental order of convergence as 

n log (E(-)/E'(.)) 
log (N h /N h ,) ’ 

where Nh and Nh/ denote the number of vertices in two consecutive meshes, with 
corresponding errors E and E'. 


5.1 Hencky-von Mises elasticity problem with analytical solution 

The first constitutive law we consider, taken from [22j, is the non-linear Hencky-von 
Mises elasticity model, for which 


er = <r(x, Vu(x)) = A(dev(£(u)))tr(e(u))/ + 2/2(dev(e(u)))e(u). 

Here above, A and p are the nonlinear Lame functions, e(u) := j(Vu + (Vu ) T ) is 
the small deformation strain tensor, the symbol tr represents the trace operator and 
dev(r) = 11 (t — ^tr(r)I) || is the Frobenius norm of the deviatoric part of the tensor r. 
We take the Lame functions as follows: 

Q Q 

p(p) := - (l + (l + /9 2 )- 1/2 )- 10 4 M P a and X(p) := - (1 - 2/r(p))-10 4 MPa Vp € M + , 

This function p corresponds to the Carreau law for viscoplastic materials. It is easy to 
verify that the hypotheses at the beginning of Section [l] are fulfilled by our choice of A 
and p. We have taken the load f such that the solution u of Problem (j2| is given by: 

u\(x, y ) = v, 2 (x, y ) = sin(7rx) sin(7ry). 


In Table [I] we report the convergence history of the virtual method (24) applied to 
our test problem with different families of meshes. The table includes the number of 
mesh vertices, the convergence rates R, and the discrete errors Eq ^ and E^ 2 - 

We observe from Table [l] that a clear first order convergence rate in the discrete H 1 
like norm and show a quadratic rate in the discrete L°° norm. 


5.2 A benchmark elasticity model problem with analytical solution 

In this test case, we select the constitutive load as 

a = tr(x, Vu(x)) = p(s( u))e(u), 
where p is defined by the following nonlinear function: 

p(e( u)) := 3(1 + || g:(u) || 2 ) • 10 4 MPa, 

with 

IH u )ll 2 = J2 Nl 2 - 

i,j =1 
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Mesh 

N h 

TT'fl 

-^0,00 

Ro,oo 

TT'fl 

^1,2 

Rl,2 


64 

3.4192e-2 

- 

4.5675e-l 

~ 


192 

8.2511e-3 

2.59 

2.4445e-l 

1.14 

njl 

640 

2.4353e-3 

2.03 

1.2803e-l 

1.07 


2304 

6.7066e-4 

2.01 

6.5274e-2 

1.05 


8704 

1.7495e-4 

2.02 

3.2919e-2 

1.03 


33792 

4.4619e-5 

2.01 

1.6527e-2 

1.02 


64 

5.6458e-2 

- 

5.0007e-l 

- 


192 

1.9675e-2 

1.92 

2.7166e-l 

1.11 


1280 

6.4750e-3 

1.85 

1.4054e-l 

1.09 


2304 

2.01403-3 

1.82 

7.1120e-2 

1.06 


8704 

5.4860e-4 

1.96 

3.5590e-2 

1.04 


33792 

1.4070e-4 

2.01 

1.7817e-2 

1.02 


25 

6.1947e-2 

- 

7.1975e-l 

- 


81 

9.3599e-3 

3.21 

3.5627e-l 

1.19 


578 

1.7576e-3 

2.62 

1.7809e-l 

1.09 


1089 

4.2329e-4 

2.14 

8.9038e-2 

1.04 


4225 

1.0516e-4 

2.05 

4.4518e-2 

1.02 


16641 

2.6254e-5 

2.02 

2.2259e-2 

1.01 


25 

1.5401e-l 

- 

1.0516e-0 

- 


81 

3.3021e-2 

2.62 

5.3972e-l 

1.14 


578 

7.1005e-3 

2.42 

2.7525e-l 

1.06 


1089 

1.6650e-3 

2.19 

1.3832e-l 

1.04 


4225 

4.1133e-4 

2.06 

6.9382e-2 

1.02 


16641 

9.0462e-5 

2.21 

3.2452e-2 

1.05 


Table 1: Approximation of u: convergence analysis of the virtual method (24). 


We have taken the load f such that the solution u of Problem Q is given by: 
ui(x,y) =U 2 (x,y) = 10sin(7rx) sin(7ry). 


We remark that this choice does not actually correspond to any elastic material. 
Instead, it has been chosen as a “benchmark model” which does not satisfy the assump¬ 
tion at the begining of Section [ij condition (28) does not hold, in particular. 

Table [2] shows the convergence history of the virtual method (24) applied to our 
test problem with different families of meshes. The table includes the number mesh 
vertices, the convergence rates R, and the discrete errors Eq^ and 2 - 

Once more, a quadratic order of convergence in the discrete L°° norm and a linear 
order convergence rate in the discrete H 1 like norm can be clearly appreciated from 
Table [2j 

We now consider the same 14 and the same constitutive law, but we choose a couple 
of different loads. The purpose is now to show the importance of updating the choice of 
the stability constant appearing in the elastic form (20), for instance by employing the 
recipe detailed in (21) (see Remark 3.2). Therefore, we consider two different external 
forces, compatible with the following two analytical solutions: 


Case 1: u = (x(l - x)y{ 1 - y),x( 1 - x)y{ 1 - y)) , 

Case 2: u = 80 * (x(l — x)y( 1 — y),x( 1 — x)y{ 1 — y)^j . 

We notice that in Case 1 the solution gives rise to deformations of moderate magnitude, 
while in Case 2 much larger deformations occur. We consider a single family of three 
regular Voronoi meshes, generated using the algorithm in [2.91. Moreover, we choose 
the following relative error measure, involving both the displacement components at all 
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Mesh 

N h 

TT'h 

-^0,00 

Ro,oo 

TT'h, 

^1,2 

Rl,2 


64 

4.1122e-2 

- 

4.6371e-0 

~ 


192 

1.7816e-2 

1.52 

2.6318e-0 

1.03 


1280 

5.0006e-3 

2.11 

1.3317e-0 

1.13 


2304 

1.2449e-3 

2.17 

6.6288e-l 

1.08 


8704 

2.9750e-4 

2.15 

3.3092e-l 

1.04 


33792 

8.2512e-5 

1.90 

1.6553e-l 

1.02 


64 

8.1685e-2 

- 

5.1698e-0 

- 


192 

2.3823e-2 

2.24 

2.9790e-0 

1.00 


1280 

1.4234e-2 

0.86 

1.5553e-0 

1.08 


2304 

5.9189e-3 

1.37 

7.6103e-l 

1.12 


8704 

1.7906e-3 

1.80 

3.6614e-l 

1.10 


33792 

4.7067e-4 

1.97 

1.7981e-l 

1.05 


25 

1.8457e-l 

- 

9.6706e-0 

- 


81 

5.2374e-2 

2.14 

4.0009e-0 

1.50 


578 

1.5787e-2 

1.89 

1.8538e-0 

1.21 


1089 

4.5978e-3 

1.86 

9.0144e-l 

1.09 


4225 

1.2340e-3 

1.94 

4.4672e-l 

1.04 


16641 

3.1086e-4 

2.01 

2.2279e-l 

1.02 


25 

1.4957e-l 

- 

11.0527e-0 

- 


81 

3.6140e-2 

2.41 

5.4418e-0 

1.20 


578 

1.1670e-2 

1.78 

2.6376e-0 

1.13 


1089 

3.6360e-3 

1.76 

1.3130e-0 

1.05 


4225 

1.1048e-3 

1.76 

6.5565e-l 

1.02 


16641 

3.1365e-4 

1.83 

3.2786e-l 

1.01 


Table 2: 


Approximation of u: convergence analysis of the virtual method (24). 


the vertices v of the mesh: 


Em = 


max ve o /t , i= i ,2 |^(v) - (u h )j(y)\ 
max ve n h , j=i ,2 K(v)| 


In Table [3] we report the relative errors computed for Case 1, using both the updated 
scalings introduced in (21) and a fixed scaling. We notice that convergence is attained 


for both the strategies of the scaling choice. 

In Table [4] we report the relative errors computed for Case 2, using both the updated 
scalings introduced in (21) and a fixed scaling. We notice that for this case, convergence 


is attained when using the updating strategy, while choosing a fixed scaling provides 
unsatisfactory results. In particular, on the finest mesh the error is still around 20%. 
Moreover, the solution is highly oscillating due to the presence of unstable numerical 
modes (figure not shown). 


Mesh 

N h 

Updated ole 

Fixed ole 

Mesh 1 

199 

1.715e—2 

1.174e—2 

Mesh 2 

800 

3.580e—3 

3.392e—3 

Mesh 3 

3179 

1.287e—3 

8.946e—4 


Table 3: Case 1: relative errors for the updated and fixed choice of the scaling. 


5.3 Von Mises plasticity 

In the present section we show a numerical example for an inelastic material, von Mises 
plasticity with linear hardenings. We consider the classical problem of a strip with 
circular hole in plain strain regime under enforced displacements of 5 amplitude at two 
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Mesh 

N h 

Updated ag 

Fixed ole 

Mesh 1 

199 

2.384e—2 

2.685e0 

Mesh 2 

800 

9.299e—3 

9.555e—1 

Mesh 3 

3179 

3.132e—3 

2.090e—1 


Table 4: Case 2: relative errors for the updated and fixed choice of the scaling. 


ends. Due to the symmetry of the problem, we can consider one quarter of the strip, 
as depicted in figure [3] (left). The geometric data are 



Figure 3: Left: depiction of the geometry for the perforated strip problem. Right: 
sample Voronoi mesh V2. 


B = 100 mm, H = 180 mm, Bq = 50 mm, 5=10 mm. 

We consider a J 2 plasticity model with linear kinematic and isotropic hardenings (see 
for instance m) with material parameters 

E = 70 MPa, v = 0.2 MPa, a y ^ = 0.8 MPa, iL; so = 10 MPa, iLkin = 10 MPa. 

For comparison purposes, we take as “exact solution” one obtained with linear finite 
elements on a fine triangular mesh with 45312 elements. Note that, since the considered 
model includes hardenings, there is no risk of volumetric locking and thus triangular 
elements are a good choice. We solve the problem on a sequence of four Voronoi meshes 
(mesh VI to mesh V4) generated with the code PolyMesher [22]. We depict a sample 
mesh V2 in figure [ 3 ] (right) while the number of vertices in each grid can be found in 


with 100 time-steps. At each time step the constitutive law is solved using a classical 
radial return map algorithm (see for instance 1223, Chapter 3). For each mesh we show 
the following values in Table [5] 


Table D} In all cases we use the incremental loading procedure described in Section 3.3 
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Mesh 

N h 

Displ. A 

Displ. B 

^max 

( Tt 

VI 

129 

0.7839 

-0.3181 

3.3842 

244.2324 

V2 

511 

0.8173 

-0.3928 

4.1354 

240.1062 

V3 

2032 

0.8253 

-0.4212 

4.4266 

238.7653 

V4 

8131 

0.8277 

-0.4300 

4.7755 

238.3688 

Reference 

22921 

0.8284 

-0.4334 

4.9891 

238.2631 


Table 5: Number of mesh vertices, displacements at points A and B, maximum stress 
and total stress for the four Voronoi meshes and for a reference value obtained with a 
fine triangular mesh. 


• The vertical displacement at the point A of coordinates (Omni, 50mm), where the 
axes origin is at the center of the hole; 

• the horizontal displacement at the point B of coordinates (50mm, 0mm); 

• the maximum stress <r max ; 

• the total stress <tt, he. the integral over O, of the stress amplitude ||cr|| = 

(Si.,-,1,2 M 2 ) 1/2 - 

Note that, on purpose, in Table [5] we consider quantities for which is easy to obtain 
convergence (displacement at point A ad total stress) and other ones for which is 
harder (displacement at point B and maximum stress). In all cases we can appreciate 
the convergence of the method towards the reference values; finer Voronoi meshes would 
be needed for a better approximation of the maximum stress. 

In figure [4] we depict the value of the plastic consistency parameter 7 for the V4 and 
for the fine reference mesh. The parameter 7 indicates if and how much plastification 
has occurred locally for the material; we refer again to {213 f° r a detailed description 
of the model. Again, the results for the proposed method are in good accordance with 
the reference one. 


5.4 Finite strain elasticity 


The method detailed in Sections |3.1||3.2 can also be applied to elastic problems in a large 
strain regime. However, we remark that the complexity of the finite elasticity problem 
requires a much deeper design and analysis than the one here presented. Therefore, 
the following discussion should be intended only as a very preliminary study towards 
the VEM discretization of large deformation elastic problems. 

We here focus on neo-Hookean hyperelastic materials, but different constitutive laws 
could be considered. Following a material description (see p [T51l25] . for instance), the 
variational formulation of the elastic large deformation problem reads as in 


Find u E V such that 


/ P(x, Vu(s)) : Vv(x) dx = [ f(x)-v(x)dx VvgW, 
JO Jfl 


(56) 


where the first Piola-Kirchhoff stress tensor P(r,Vu(x)) is not necessarily symmetric. 
As for Problem ([ 3 ]), in (56) the symbol V denotes the space of admissible displacements 
and W the space of its variations. A homogeneous neo-Hookean material is described 
by the constitutive law: 
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Figure 4: Depiction of the plastic flow 7 , mesh V4 on the left and reference triangular 
mesh on the right. 


P{x, Vu(x)) =m[(I + Vu) + (I + Vu)" T ] 

+ A0(det(I + Vu)) 7 r(det(I + Vu))(/ + Vu)“ T . 


(57) 


Above, A and ^ are given constants, 0 : M + —M 
7 r is defined as 

7r(s) = 0 / (s)s . 


is a suitable smooth function, and 

(58) 


Here, we choose 0(s) = s — 1, so that 7 r(s) = 1. 

A possible virtual method for Problem (56) can be designed exactly as in Sec¬ 


tions 3.1 3.2 , simply by systematically substituting P in place of er. 


We test the method considering a square block of side length lrn, which initially 
occupies the region D = (0, l) 2 . We impose clamped boundary conditions on the side 
r c = { 0 } x [0,1], while the remaining part of the boundary is free. The material 
parameters are chosen as ^ = 2.6316 • 10 4 MPa and A = 5.1086 • 10 4 MPa. The load is 
given by f = (1,0) T 10.5 • 10 10 N/m 3 . 

Table [b] displays the computed displacements of the material point P = (1, 1) T , 
when using triangular (T1,...,T4), quadrilateral (Q1,...,Q4), and hexagonal Voronoi 
(V1,...,V4) meshes. A reference solution at the same point, obtained with a very fine 
triangular mesh of 70344 elements, corresponding to 35459 mesh vertices, is also re¬ 
ported. Finally, Figure [5] depicts the deformed body when using the triangular mesh 
T2, the square mesh Q2 and the hexagonal Voronoi mesh V2 of Table [ 6 j We notice 
that for every considered scheme, convergence to the reference solution occurs, and the 
deformed shapes appear to be sensible. 


6 Conclusions 

We have presented a Virtual Element Method to deal with fairly general non-linear 
elastic and inelastic problems. Our scheme is based on a low-order approximation of the 
displacement field, together with a suitable treatment of the numerical displacement 
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Mesh 

N h 

x— Displ. at P 

y— Displ. at P 

T1 

55 

0.9865 

-0.0438 

T2 

183 

1.0615 

-0.0398 

T3 

727 

1.0848 

-0.0358 

T4 

2810 

1.0967 

-0.0354 

Q1 

49 

0.9979 

-0.0736 

Q2 

196 

1.0730 

-0.04791 

Q3 

784 

1.0950 

-0.0391 

Q4 

3025 

1.1005 

-0.0364 

VI 

52 

0.9125 

-0.0673 

V2 

199 

1.0344 

-0.0520 

V3 

800 

1.0722 

-0.0408 

V4 

3179 

1.0918 

-0.0368 

Reference 

35459 

1.1018 

-0.0353 


Table 6: Computed displacements using triangular (T1,...,T4), square (Q1,...,Q4), and 
hexagonal Voronoi (V1,...,V4) meshes. 



Figure 5: Deformed body obtained with the triangular mesh T2 (left), the square mesh 
Q2 (center), and the hexagonal Voronoi mesh V2 (right). 


gradient. The proposed method allows for general polygonal/polyhedral meshes, is 
efficient in terms of number of applications of the constitutive law, and can make 
use of any standard black-box constitutive law algorithm. We have presented several 
numerical tests assessing the computational performance of the proposed methodology. 
However, we remark that this study is intended as a first step towards the design of 
efficient Virtual Element Methods for non-linear Computational Mechanics problems. 
Many possible extensions and improvements could be of interest. For instance, large 
deformation problems require a much deeper investigations, and other inelastic cases 
such as perfect plasticity or damage could be considered. 


Acknowledgements. D. Mora was partially supported by CONICYT-Chile through 
FONDECYT project No. 1140791 and by project Anillo ACT 1118 (ANANUM). 
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